Journal of Theoretical Biology

Global populations are ageing at an unprecedented rate. For many diseases, age is a strong indicator of susceptibility, morbidity, or mortality. Principles of evolutionary biology can be leveraged to understand how pathogens may optimally exploit new populations, and the impact of this on the global burden of infectious-disease-induced mortality. We parameterised an age-specific R0 model with 2017 epidemiological data on Measles, Tuberculosis, Meningitis, and Ebola, and age-specific demographic estimates for 2017 and 2050, for the seven Global Burden of Disease super-regions. We explored the theoretical trade-offs between pathogen virulence & transmission, and virulence & host recovery, parameterising trade-off parameters using Latin Hypercube Sampling. Population ageing between 2017-2050 saw an increase in virulence induced mortality in four settings: 1) Ebola in sub-Saharan Africa, 2) Measles in central/eastern Europe & central Asia region, 3) Measles in North Africa & the Middle East and 4) Tuberculosis in the central/eastern Europe & central Asia region. The decrease in infection duration due to an increase of elderly people drives pathogen virulence down for diseases in the remaining settings. Understanding the mechanisms that shape pathogen dynamics and leveraging this to predict future challenges is key to endemic disease management in a rapidly changing world. Author SummaryKey aspects of disease transmission including susceptibility to infection, the severity of infection, and the probability of dying from that infection, are affected by host age. Global populations are rapidly ageing, so that the mean age of most populations is generally higher than it used to be and is set to continue on this trajectory. This suggests that the dynamics of infectious diseases are also likely to change, although infectious disease dynamics tend to be non-linear as these key parameters interact. We have developed a dynamic modelling framework to explore how changes in population age structure might impact the optimal level of pathogen virulence in a population. We have chosen four infectious diseases as case studies, that differentially impact certain age classes to illustrate these dynamics. We have parameterised this framework with open access data for each of the seven Global Burden of Disease super-regions and show that population ageing can increase virulence for several diseases in differing global regions, whilst increased background rates of mortality can drive virulence down in others.

Trees accumulate somatic mutations throughout their long lifespan, resulting in genetic mosaicism among branches. While recent genomic studies quantified these mutations, they were largely limited to describing static patterns of variation. In this study, we developed a mathematical model to infer the dynamic processes of somatic mutation accumulation from snapshot genomic data obtained from four tropical trees (Dipterocarpaceae), which dominate tropical rain forests in Southeast Asia. Our model focus on genetic differences between shoot apical meristems (SAMs) at branch tips and explicitly incorporate stem cell dynamics within SAMs during shoot elongation and branching, enabling us to quantify somatic genetic drift arising from stem cell lineage replacement. By comparing model predictions with empirical data from Dipterocarpaceae trees, we estimated key parameters governing stem cell dynamics and somatic mutation rates. Our results indicate that both shoot elongation and branching involve replacement of stem cell lineages, leading to a moderate degree of somatic genetic drift. Accounting for stem cell dynamics resulted in slightly lower mutation rate estimates than previous approaches that ignored these processes. Using the estimated parameters, we further performed stochastic simulations to predict patterns of somatic mutations, including features not directly observed in the sampled trees, such as occasional deviations of somatic mutation phylogenies from physical architecture. Together, our modeling framework provides insights into how genetic mosaicism is shaped within tropical trees and reveals the stem cell dynamics underlying their long-term growth and accumulation of somatic mutations. (236 words) Highlights- We built mathematical models to predict the genetic differences between branch tips by somatic mutations. - The model considers the varying dynamics of stem cells in shoot meristem during shoot elongation and branching. - We compared the model prediction with empirical data from tropical trees, Dipterocarpaceae, and estimated the dynamics of stem cells and mutation rate. - Somatic mutation dynamics were shaped by somatic genetic drift arising from stem cell lineage replacement during shoot elongation and branching. - Accounting for stem cell dynamics led to slightly smaller estimates of mutation rates compared with previous estimates that ignored the dynamics. - Our models offer insights into how genetic variability is shaped in the tropical trees and the stem cell dynamics underlying their long-term growth.

Alzheimers disease (AD) is characterized by the accumulation of amyloid-{beta} (A{beta}), yet the specific link between plaque burden and cognitive decline remains a subject of intense investigation. This paper presents a mathematical model that simulates the coupled dynamics of A{beta} monomers, soluble oligomers, and fibrillar species in the brain tissue. By modifying existing moment equations to include a dedicated conservation equation for A{beta} monomers, the model explores how various microscopic processes, such as primary nucleation, surface-catalyzed secondary nucleation, fibril elongation, and fragmentation, contribute to macroscopic disease progression. Central to this study is the concept of "accumulated neurotoxicity" as a surrogate marker of biological age, defined as the time-integrated concentration of soluble A{beta} oligomers. Unlike plaque burden, accumulated neurotoxicity cannot be reversed, and the harm it causes depends critically on the sequence of events that produced it. Numerical results demonstrate that while plaque burden and neurotoxicity both increase over time, their relationship is non-linear and highly sensitive to the efficiency of protein degradation machinery. Specifically, impaired degradation leads to a rapid advancement of biological age relative to calendar age. The model further identifies oligomer dissociation and fibril fragmentation as potential protective mechanisms that can counterintuitively reduce neurotoxic burden by diverting monomers away from the soluble oligomer pool. These findings provide a quantitative framework for understanding why individuals with similar plaque burdens may experience vastly different cognitive outcomes, underscoring the importance of targeting soluble oligomers early in therapeutic interventions.

In this research, we create a new fractional-order SEIHRD framework to examine how the Nipah virus moves from one species to another (zoonotic spillover) and how it later spreads throughout a community (via contact with one another) or in a hospital or isolation situation (via entering into a hospital or being placed under quarantine). We used the fractional-derivative formulation of the SEIHRD model to demonstrate memory-based effects related to the progression of an infection and also reflect time-distributed effects associated with surveillance and control measures placed on an infected patient. We first demonstrated that the basic epidemiologic properties of the model were consistent by showing that the solutions of the SEIHRD differential equations will always yield positive and bounded solutions within biologically relevant parameter ranges. We then established the well-posedness of this model by transforming the SEIHRD differential equations into an equivalent integral operator and applying various fixed-point arguments to demonstrate that there will always be unique solution(s) to the SEIHRD differential equations. To evaluate the threshold parameter for the transmission of Nipah virus within a given population we calculated the threshold level through the next generation method to determine the expected number of secondary infections from a new or chronically infected host. One of the main contributions of this work is to include an analysis of the robustness of a given solution to all potential perturbations (i.e., Ulam-Hyers and generalized Ulam-Hyers stability). In addition, we provide analytic results guaranteeing that small perturbations due to approximate modeling, numerical approximation (discretization), or the lack of data fidelity will produce controlled deviations in the solutions. To finish this project, we perform a global sensitivity analysis on uncertain coefficients to evaluate their contribution to the uncertainty of each coefficient and to find out the coefficients that most strongly influence major outcome metrics. This will allow us to develop a priority order for prioritizing spillover control (reduction of human contact and/or isolation), contact reduction, and expenditure of resources towards isolation-related interventions. The resulting framework converts fractional epidemic modeling from a descriptive simulation to a replicable method with robustly defined behavior and equal response prediction.

Developing effective antimicrobial strategies requires a predictive understanding of bacterial responses to multiple stress conditions which often result in multiple phenotypes. A microbes survival and proliferation depend on its ability to manage concurrent, dynamically varying stressors within its microenvironment. However, time-resolved predictive models that capture multi-phenotype responses are lacking, and single-phenotype models often fail to accurately replicate a microbes reaction to mixed stress conditions. In this work, we develop a mechanistic in-silico model of multi-stress response in Escherichia coli K12 and use it to characterize phenotype dynamics in iron-limited and hydrogen peroxide containing environments. Specifically, we replicate the iron and oxidative stress response networks in E. coli using a system of ordinary differential equations and applied a multi-phenotype parameterization scheme that leverages multi-measure empirical data, augmented metric-based sensitivity analysis, sequential parameter estimation, and ensemble modeling. Our approach resulted in robust models with a 93% accuracy when compared to experimental datasets across 20 stress-response categories, outperforming traditional single-phenotype approaches (80-87%). Analysis of posterior parameter distributions revealed that multi-phenotype optimization eliminates heavy-tailed distributions characteristic of poorly constrained fits and shifts parameter posteriors from boundary-concentrated to centrally localized forms, indicating improved identifiability. Simulation outcomes confirmed key features of E. colis iron metabolism, showing that moderate peroxide stress in an iron-rich environment creates significant adaptation challenges, leading to a bacteriostatic phenotype. The model provides insights into biochemical mechanisms important to E. colis temporal response to varying iron availability, with implications for ecological dynamics and pathogenesis. Our parameterization approach highlights the effectiveness of a combination of optimization methods and ensemble modeling in developing predictive models that are robust across multiple phenotypes. Results demonstrate that data structure, specifically the integration of multiple phenotypes and response outputs, proves to be as critical if not more critical than data volume for achieving well-constrained parameter estimates and robust predictions across experimental conditions. Author SummaryUnderstanding bacterial stress response is crucial for developing strategies to control bacterial populations, particularly as antibiotic resistance poses a growing threat, with "superbugs" potentially triggering a global health crisis. While mathematical models offer powerful tools to study biological systems, many struggle to predict cellular behavior across multiple phenotypes due to the complexity of responses. Iron metabolism is vital for bacterial survival, particularly under oxidative stress, leading to various bacterial growth dynamics. This work uses mathematical modeling to explore how E. coli manages multiple stressors, focusing on iron metabolism and oxidative stress. By applying a novel combination of optimization and ensemble modeling methods, we improved model accuracy by nearly 16%, enabling predictions of E. colis varied response to single, dual, and dynamic stress environments. Our approach offers a valuable tool for understanding and combating bacterial persistence, with future studies able to expand its use to determine how bacterial communities respond to multiple stressors.

The emergence of antimicrobial resistant bacteria has been identified as one of the most serious public health and development threats for the near future. The use of bacteriophages (phages) is a promising solution for the sustainable control of these pathogens. Phages are natural viral predators of bacterial pathogens. However, due to the variability and adaptability of bacteria, developing effective and sustainable phage treatments requires drawing from a wide variety of different phage species. This study applies specialised mathematical and computational estimation approaches to the problem of sampling and discovering species of phages in microbiological communities. We show that classical non-parametric estimator techniques lead to robust results and outperformed, for existing data settings in phages, model-based approaches. We then show how efficient the continuation of current phage collection and isolation effort is expected to be for discovering new phage species in various relevant bacterial host genera, a prerequisite for phage applications to provide sustainable control of pathogens in human and environmental settings. Our results have the potential to inform and optimise the hunt for and isolation of novel of phages from the the natural environment.

Mammalian spermatogenesis occurs in the seminiferous tubules, which exhibit unique spatiotemporal differentiation patterns known as cellular association patterns. In mice, these patterns can be regarded as one-dimensional wavetrains that consistently propagate inward from both ends, resulting in one or more "sites of reversal." Segmented wavetrain pattern, in which the wave propagation direction spatially switches, was observed in our previous three-species reaction-diffusion model for interspecific species difference in spermatogenic waves (Kawamura et al., 2021). However, the biological mechanisms of the formation of sites of reversal and of this directional bias, as well as the principle of pattern formation, remain unknown. Here, we refined our previous model to match the actual biological spatiotemporal scale and examined its dynamics through extensive numerical simulations. The modified model frequently generated segmented wavetrain patterns, corresponding to the sites of reversal, but without directional bias. We systematically examined possible biological mechanisms for the bias and found that tubule elongation, especially near the rete testis, most effectively accounts for the bias among the tested. Extensive simulations revealed that the segmented pattern is numerically stable, emerges more frequently in longer domains, and shows an exponential segment size distribution with a lower limit for the stably existing segment length. These explorations imply that locally emerged unidirectional wavetrains serve as building blocks to generate the stable segmented wavetrains through their interactions. HighlightsO_LISegmented wavetrains reflect sites of reversal in seminiferous tubules. C_LIO_LISegmented patterns frequently emerge but show no inherent directional bias. C_LIO_LITubule elongation may contribute to inward propagation near the rete testis. C_LIO_LISegmented wavetrains are numerically stable and more frequent in longer domains. C_LIO_LIInteractions of local unidirectional wavetrains generate stable segmented structures. C_LI

Studies of appendage regeneration in vertebrates have shown that the fundamental building block of any regenerative tissue is a blastema. The blastema is a cone-shaped accumulation that forms at the site of amputation post wound-healing and is the result of a highly coordinated process involving a cluster of cells capable of growth, migration, and differentiation. Although several key signaling pathways involved in regeneration have been identified, which cellular processes they control and how these processes are coordinated across space and time are not yet fully understood. This study introduces a computational tool to examine how the outgrowth results from the interaction of two tissue layers: the bulk (mesenchyme) and the overlying epithelium. We developed a novel hybrid agent-based modeling framework and an accompanying parameter inference pipeline to uncover the cellular properties in the epithelium and the mesenchyme driving the formation of a normal regenerative blastema with a morphology similar to that observed in experiments. Using our model, we report two conditions for blastema formation: retained local softening of the epithelial layer at the site of injury, which was confirmed experimentally with atomic force microscopy (AFM) measurements, and the involvement of the Wnt signaling pathway in the directed migration of mesenchyme cells towards the distal tip. Taken together, this combined experimental-theoretical approach provides a framework for understanding how the Wnt signaling pathway influences the formation of the early blastema at multiple levels of organization and how key cellular behaviors contribute to its formation. Author SummaryA small number of tetrapods have retained the ancestral ability to regenerate tissues and even limbs. Indifferent of species or tissue, the decisive initial stage of limb regeneration is the formation of a specialized structure called the blastema, a heterogeneous mass of mesenchymal cells, in a relatively short timescale of 2-14 days post injury or amputation. To study the mechanical and cellular conditions for limb blastema formation in the axolotl model organism, we developed a novel hybrid agent-based modeling framework and accompanying kinetic parameter inference pipeline. By recapitulating blastema morphometrics of healthy and stalled regenerative states, our model finds two conditions for blastema formation: retained local softening of the epithelial layer at the injury site post wound-healing, which we confirmed with atomic force microscopy measurements, and that the Wnt signaling pathway plays a role in the migration of mesenchyme cells to the distal tip in order to produce the blastema.

Evolutionary cancer therapy (ECT) applies principles of evolutionary game theory to prolong the effectiveness of cancer treatment by curbing the development of treatment resistance. It was shown to increase time to progression while decreasing the cumulative drug dose. ECT individually tailors treatment schedules for patients based on their cancer dynamics and, thus, requires regular follow-up and precise measurements of the cancer burden. The current literature on ECT often overlooks clinical realities, such as rather long intervals between tests, possible appointment delays and measurement errors, in the development of the treatment protocols. In this study, we assess the clinical feasibility of ECT for metastatic non-small cell lung cancer (NSCLC). We create virtual patients with cancer dynamics described by the polymorphic Gompertzian model, based on data from the START-TKI clinical trial. We assess the effects of longer test intervals, measurement error and appointment delays on the expected time to progression under the evolutionary therapy protocols. We show that a higher containment level, although it increases time to progression in the models predictions, may lead to premature treatment failure in the presence of measurement error and appointment delay. Further, we show that the ECT protocol with a single containment bound is more robust to the clinical realities than the protocol with two bounds. Finally, we show that a dynamically adjusted treatment protocol can be beneficial for individual patients, but requires a thorough follow-up. This study contributes to the design of a clinical trial and the future clinical implementation of evolutionary therapy for NSCLC.

Patterns in the vegetation across arid and semiarid regions may be explained as a form of self-organization driven by water scarcity, and are often modeled through reaction-diffusion dynamics. Recent work has shown that similar mathematical models generate patterns on networks. However, these studies have focused on idealized topologies with no reference to natural pattern-forming systems. Our study aims at bridging these two fields: we employ a physical reaction-diffusion vegetation model, and gradually modify the topology of the diffusion network by adding random shortcuts over a 2-dimensional grid, interpolating between a regular lattice and a random network. We found that network topology strongly shapes both the resulting vegetation patterns and the precipitation range that supports them. Three behavioral regimes emerge. On a regular lattice, high-regularity patterns develop reflecting local diffusion processes. On a random network, the system is dominated by global pressure towards homogenization yielding either a uniform state or a single patch. In the intermediate shortcut density range, as the network topology resembles a small world network, the interaction between the two scales of diffusion generates two kinds of disordered patterns: low-regularity patterns with a well-defined characteristic wavelength, and irregular patterns characterized by a broad patch size distribution. These disordered patterns resemble real-world observations and, in our model, they show different responses to changing precipitation. Although we focused on dryland vegetation, we suggest that network-mediated diffusion could lead to similar mechanisms in a wide variety of pattern-forming systems. HighlightsO_LIWe study vegetation pattern formation over different diffusion network topologies. C_LIO_LITwo kinds of stable disordered patterns states develop over small world topologies. C_LIO_LILow-regularity patterns with a well-defined characteristic wavelength. C_LIO_LIIrregular patterns characterized by a broad patch size distribution. C_LIO_LIThese different kinds of disordered states show different relations to precipitation. C_LI

Given the accumulation of evidence that evolution can affect ecological dynamics, especially under global change scenarios, a key question is how such ecoevolutionary dynamics may change the coexistence of species and biodiversity in general. In the present article, we propose a graphical approach allowing to simultaneously discuss ecological coexistence and phenotype evolution. Our graphical approach allows tackling the two aspects in the same parameter space, allowing direct links between ecological and evolutionary perspectives. While evolution is often thought positive for the resilience of ecological systems, we first highlight it does not usually allow for better coexistence for the system as a whole. Even when focusing on the fate of the species that is evolving, evolution often leads to greater vulnerability. The graphical approach we propose is flexible and can be applied to all interaction types and covers variations in trade-off structures. Using this flexibility, we highlight how evolutionary effects can be positive or negative for coexistence, depending on these two components. Finally, we illustrate how the approach can be applied, using empirical examples derived from the literature. We thereby highlight the critical ingredients needed to inform the graphical approach, its potential use for proposing testable scenarios, but also clarify its limits.

Our recent work on molecular evolution and population genetics postulated that individuals with a specific mutation exhibit a fluctuation in fitness, short for FSI (fluctuating selection among individuals), whereas the fitness effect of wildtype remains a constant. An intriguing phenomenon called selection-duality emerges, that is, a slightly beneficial mutation could be a negative selection (the substitution rate less than the mutation rate). It appears that selection-duality is bounded by two bounds: the generic neutrality where the mutation is neutral by the means of fitness on average, and the substitution neutrality where the substitution rate equals to the mutation rate. In addition, the middle point of generic neutrality and substitution neutrality is called the FSI-neutrality. An important problem is about the age profile of allele frequency, i.e., the arising timing of a mutation whose frequency in the current population is given (the allele-age problem for short). Solving this problem under selection duality would help extend the standard coalescent theory that based on strict neutrality to a more general form under selection duality. In this paper, we studied the allele-age problem under selection-duality by the first arrival time approach and the mean age approach, respectively. Since the general solution of allele-age problem under selection duality is not available, we focused on solving the problem at the substitution neutrality (the up-bound of selection duality), the FSI-neutrality (the middle-point) and the generic neutrality (the low-bound), respectively. Our analysis results in an overall picture that the mean first-arrival age of a mutation at the substitution neutrality is theoretically identical to that at the FSI-neutrality, which is numerically close to that at the generic neutrality. For illustration, we calculated the mean age of nonsynonymous mutations in the human population and demonstrated that the estimated allele-age could be overestimated considerably when the effect of FSI was neglected.

The transbilayer distribution of plasma membrane cholesterol remains uncertain despite repeated analysis. We propose a new mechanism driving cholesterol sidedness: sterols form simple stoichiometric associations with phospholipids. Our model postulates that the phospholipids in the plasma membrane bilayer are fully complexed with cholesterol. The cholesterol in each leaflet is then the product of the abundance of its phospholipid and its sterol stoichiometry. Notably, lipid affinities are not relevant. Applying literature values for the composition, abundance and sterol stoichiometry of the phospholipid in each leaflet, the model predicts that two-thirds of the cholesterol in the human erythrocyte membrane bilayer is located in its outer leaflet, an exofacial to endofacial ratio of 2:1. The model also predicts that the overall cholesterol content of the bilayer is [~]0.75 mole/mole phospholipid, in agreement with literature values. Furthermore, our analysis suggests that the areas of the two membrane leaflets are about the same. The concordance of prediction with observation validates the model and the values used for the parameters. The sterol in the exofacial leaflet of the plasma membrane of any cell is predicted to exceed that on its contralateral side when its phospholipids have a higher sterol stoichiometry and are fully complexed. SynopsisWe propose that the transbilayer distribution of cholesterol in the plasma membrane bilayer is determined by its complexation with the phospholipids in the two leaflets. Because the complexes are homeostatically filled to stoichiometric equivalence, leaflet cholesterol is given by the abundance of its phospholipids multiplied by its sterol stoichiometry. The model predicts that two-thirds of the cholesterol in the human erythrocyte membrane bilayer resides in the outer leaflet. It also predicts the cholesterol content of the bilayer as a whole.

Insulin-stimulated glucose uptake (ISGU) in adipocytes is central to maintain systemic glucose homeostasis. Understanding how ISGU relates to adipocyte traits, such as cell-size, is critical for elucidating pressing questions related to metabolic dysfunction connected to adipose hypertrophy and hyperplasia. Cell size is considered a central trait reflecting multiple aspects of adipocyte metabolism. However, a robust quantitative approach to estimate ISGU for a specific cell size is currently missing. Here, in an attempt to move towards such a method, we have formulated an approach using a mathematical framework. The framework consists of a linear equation: the product of the known number of cells (calculated using coulter counter-based cell-size distributions) and the unknown ISGU/cell, compared to the absolute ISGU (measured using 14C-glucose tracer assays). To solve this equation, we formulate a minimization problem which is optimized to find the unknown ISGU/cell for the best solution. Using different formulations of the equation we can investigate the need for either cell size-dependent or independent ISGU/cell, to describe varying glucose uptake in a cell sample of various cell sizes. While the framework needs further refinement, we demonstrate that cell size-dependent uptake slightly improved the agreement between model and experimental data for some groups. Together, with further validation this could serve as a useful tool to resolve long-standing questions regarding size-dependent characteristics like adipocyte size and cellular function. Key findingsHerein we explore a method to quantify cell size-dependent glucose uptake in adipocytes

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

The BCG vaccine remains the only licensed vaccine against tuberculosis (TB), yet the mechanisms behind BCG-induced protection remain poorly understood. Plumlee et al. (PLOS Pathogens 2023) infected over 1,000 mice, half of which were vaccinated with BCG, with an ultra-low dose (ULD) of Mycobacterium tuberculosis (Mtb); the authors found that BCG vaccination resulted in fewer infected mice, lower CFU lung burden, and more frequent unilateral lung infection. We have developed several mathematical models of Mtb dynamics and dissemination between murine right and left lungs and fit these models to the CFU data from unvaccinated or BCG-vaccinated mice. Alternative mathematical models incorporating either direct (lung-to-lung) or indirect (lung-intermediate-tissue-lung) dissemination pathways fit the unvaccinated data equally well, suggesting multiple plausible routes of Mtb spread. Yet, irrespective of the dissemination route, the models predicted rapid Mtb replication during early infection, transient control within 1-2 months after infection, and continued bacterial growth in the chronic phase. Fitting models to the data from BCG-vaccinated animals revealed that BCG reduces the rate of Mtb dissemination between the lungs by 89% while having a more modest effect on the replication rate within the lung, reducing it by 9%. We found that the dominant effect of BCG in curbing lung dissemination arises from its ability to reduce Mtb replication resulting in fewer infected mice, lower lung CFU, and decreased bilateral infection of the lung. We used our parameterized mathematical models to calculate the number of mice needed to detect the efficacy of a hypothetical vaccine on the probability of Mtb clearance or dissemination between murine lungs that extends previously provided estimates. Taken together, our novel mathematical modeling-based framework provides a rigorous way of quantifying vaccine efficacy in ULD-infected mice, paving the way for the pre-clinical evaluation of next-generation TB vaccines. Author summaryThe Bacillus Calmette-Guerin (BCG) vaccine remains the only licensed vaccine against tuberculosis (TB), a disease caused by Mycobacterium tuberculosis (Mtb) bacteria. BCG is clearly protective for several years when given at birth but the mechanisms by which it provides protection remain incompletely understood. By combining analysis of Mtb dynamics in over 1,000 control and BCG-vaccinated mice, infected with an ultra low dose of Mtb, with mathematical models we show that BCG reduces both the rate of Mtb replication in the lung and the rate of Mtb dissemination between the lungs. Importantly, our model suggests that BCG may thus reduce a chance of being infected and that BCG-mediated prevention of disseminated TB in humans could be due vaccines ability to block Mtb dissemination out of the lung. Our novel framework thus will allow to predict efficacy of the next generation TB vaccines in settings of infection at doses that humans are thought to be typically exposed to.

Mass-conserved reaction-diffusion systems are used as mathematical models for various phenomena such as cell polarity. Numerical simulations of this system present transient dynamics in which multiple stripe patterns converge to spatially monotonic patterns. Previous studies indicated that the transient dynamics are driven by a mass conservation law and by variations in the amount of substance contained in each pattern, which we refer to as "pattern flux". However, it is challenging to mathematically investigate these pattern dynamics. In this study, we introduce a reaction-diffusion compartment model to investigate the pattern dynamics in view of the conservation law and the pattern flux. This model is defined on multiple intervals (compartments), and diffusive couplings are imposed on each boundary of the compartments. Corresponding to the transient dynamics in the original system, we consider the dynamics around stripe patterns in the compartment model. We derive ordinary differential equations describing the pattern dynamics of the compartment model and analyze the existence and stability of equilibria for the reduced ODE with respect to the boundary parameters. For a specific parameter setting, we obtained results consistent with previous studies. Moreover, we present that the stripe patterns in the compartment model are potentially stabilized by changing the parameter, which is not observed in the original system. We expect that the methodology developed in this paper is extendable to various directions, such as membrane-induced pattern control.

Human walking is intrinsically variable. For example, there is considerable stride to stride variability even when walking speed is constant. This variability is due to uncertainty in the sensorimotor system and the environment, and is shaped by both musculoskeletal dynamics (e.g. joint stiffness and damping originating from muscles) and the control strategy used to mitigate the effects of uncertainty. Yet, insight into how sensorimotor noise shapes walking variability is limited due to a lack of experimental methods to assess sensorimotor noise and control strategies during walking. Simulations that account for uncertainty can elucidate how sensorimotor noise affects movement variability but due to numerical challenges, accounting for sensorimotor noise is not common in simulations of walking. Existing simulations have hugely simplified musculoskeletal dynamics (e.g. no muscles), the control policy (e.g. pre-defined feedback loops), or sensorimotor noise sources (e.g. only motor noise). Here, we performed stochastic optimal control simulations of walking based on a model with 9 degrees of freedom and 18 muscles to study how the level of sensory and motor noise influences walking. We solved for feedforward muscle excitations and full-state time-varying feedback gains that minimised expected effort while generating periodic, and hence stable, gait patterns. To enable these simulations, we approximated the state distribution with a Gaussian and used an unscented transform to propagate the state covariance. Resulting optimisation problems were solved with direct collocation. Sensorimotor noise level had a small effect on the mean kinematics but shaped kinematic and muscle activity variability as well as expected effort. Although simulations underestimated the magnitude of experimental positional variability, they captured its structure. In agreement with experimental results, the control policy prioritised limiting variability of centre of mass kinematics and minimal swing foot clearance over limiting joint angle variability. Hence, our simulations suggest that effort minimisation underlies these observations. Author summaryWhen performing a movement multiple times, each repetition will be slightly different due to random disturbances in the neural signals used to control movement, i.e. sensorimotor noise. Because it is difficult to measure inside the nervous system of a moving person, computer simulations are used to study movement control. They found that both sensorimotor noise and musculoskeletal mechanics determine how people control arm movements and standing. However, there are no simulations of walking that systematically evaluated how sensorimotor noise level influences walking kinematics because they pose computational challenges. Here, we proposed and used an approach for minimal effort simulations of walking in the presence of uncertainty. We imposed forward speed and stability but not kinematics. We found that the level of sensorimotor noise had little effect on the mean movement but a strong effect on the variability and the expected effort. The control strategy prioritised reducing the variability of the centre of mass position and swing foot clearance over reducing the variability of individual joint angles, which is also observed in experiments. Interestingly, strict control of centre of mass position and foot clearance in our simulations emerged from minimising effort.

Gene tree parsimony (GTP) is a common approach for efficient reconciliation of multiple discordant gene tree phylogenies for the inference of a single species tree. However, despite the popularity of GTP methods due to their low computational costs, prior work has shown that some commonly employed parsimony costs are statistically inconsistent under the multispecies coalescent process. Furthermore, a fine-grained analysis of the inconsistency has indicated potentially complimentary behavior of duplication and deep coalescence costs for symmetric and asymmetric species trees. In this work, we prove inconsistency of GTP estimators for all linear combinations of duplication, loss and deep coalescence scores. We also explore empirical implications of this result evaluating inference results of several GTP cost schemes under varying levels of incomplete lineage sorting.

AA amyloidosis is a severe complication of chronic inflammatory diseases characterized by fibrillar protein deposition in the kidneys, leading to progressive organ failure. This study presents a mathematical model coupling SAA-HDL binding dynamics with renal amyloid aggregation kinetics to elucidate disease pathogenesis. Under normal conditions, Serum Amyloid A (SAA) circulates bound to high-density lipoprotein (HDL), which acts as a molecular chaperone preventing misfolding. However, during chronic inflammation, SAA production exceeds HDL binding capacity, resulting in free SAA that undergoes renal filtration. The model calculates free SAA concentration from reversible binding equilibrium and incorporates renal filtration, mesangial accumulation, and conversion to amyloid fibrils through primary nucleation and autocatalytic growth mechanisms. A central contribution of this work is quantifying accumulated nephrotoxicity arising from AA oligomers, which inflict cumulative cytotoxic damage to mesangial and tubular cells over time. Because oligomers are continuously generated during ongoing aggregation, their toxic burden integrates across the entire duration of the disease. Combined nephrotoxicity, encompassing both oligomer-mediated cellular injury and fibril-driven mechanical disruption of renal architecture, therefore reflects not merely the current disease state but the full inflammatory trajectory of the patient. This cumulative damage defines renal biological age, a measure of functional deterioration whose portion attributable to accumulated nephrotoxicity is irreversible. Renal biological age is also path-dependent: two patients with identical present-day SAA levels may carry different renal damage burdens depending on the duration, timing, and severity of their prior inflammatory episodes. Sensitivity analysis reveals that HDL concentration and SAA cleavage rate are critical determinants of amyloid burden.